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t-h Abstract 

We present a novel method to study the dynamics of bulk fermion systems such as the neutron- 
star crust. By introducing periodic boundary conditions into Fermionic Molecular Dynamics, it 
becomes possible to examine the long-range many-body correlations induced by antisymmetrisa- 
tion in bulk fermion systems. The presented technique treats the spins and the fermionic nature 
of the nucleons explicitly and permits investigating the dynamics of the system. Despite the in- 
creased complexity related to the periodic boundary conditions, the proposed formalism remains 
computationally feasible. 
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1 Introduction 

^ J Among the myriad of phenomena in the universe, a neutron star is one of the extremely interesting ob- 

jects. As an "oversized nucleus", it is one of the most dense and compact objects known. The densities 

^ °f the inner part of a neutron star are expected to be as large as five to ten times the normal nuclear 
saturation density n . This leads to intriguing speculations about the composition of the neutron-star 
core pQ. The densities of the outer layers of a neutron star are on the average below n . At these 
densities, the neutron-star crust consists of protons, neutrons and electrons, organising themselves as a 

03 result of the short-range nuclear attraction and the long-range Coulomb repulsion. 



In the outer crust of a neutron star, i.e. the lower-density regions, nucleons bind in nuclei and are 
placed on a Coulomb lattice, embedded in an electron gas. These nuclei grow with increasing density. 
When entering the inner crust, the nuclei surpass the neutron-drip densities and immerse the Coulomb 
lattice in a neutron gas. Finally, in the densest layers of the crust (no/3 to n /2 [lj), the spherical 
nuclear clusters outgrow the unit cell of the Coulomb lattice, transforming the crystallised crust into 
uniform core matter. In this crust-core interface, also known as the mantel of the neutron star, the sys- 
tem balances on a subtle interplay between the nuclear surface energies and the Coulomb energy of the 
neutron-proton-electron system. This results in a multitude of competing quasi-ground states, having 
sometimes quite different matter distributions, where the system can cool down to. These complex- 
shaped structures are dubbed "nuclear pastas" [2J[3llI]. 
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The crust-core interface of a neutron star may represent a sizable fraction of the crustal mass [5]. 
After all, it represents the densest regions of the crust. Hence, to apprehend the crust is an important 
key to understanding many astrophysical phenomena related to neutron stars. Furthermore, matter at 
sub-saturation densities, including nuclear pastas, can contribute 10-20% of the mass in the later stages 
of a collapsing stellar core [6j, making them also important for supernovae physics. As such, a sound 
description of the neutron-star crust, or nuclear matter at sub-saturation densities in general, is needed 
to understand various astrophysical observations [4J. 

The nuclear pasta phases have been investigated with a broad spectrum of computational tech- 
niques ranging from the liquid drop model over Hartree Fock to Monte Carlo and molecular dynamics 
El El Dm CD]. From all these, molecular dynamics is one of the few that is not biased with regard 
to the geometry of the nuclear clusters. Moreover, thanks to the dynamical aspect, one can go beyond 
the study of ground-state structures. The latter is most interesting for nuclear pastas as, with a pre- 
ponderance of low-energy excitation levels, they are susceptible to low-energy dynamics stemming from 
external probes or temperature changes. 



2 Fermionic Molecular Dynamics 

Molecular dynamics techniques have already been successfully deployed to study various properties of 
nuclear pastas and the neutron-star crust. Classical Molecular Dynamics (CMD) provided a way to 
study the effect of large density fluctuations on neutrino opacities and to evaluate the breaking strain 
of the crust |T2l [13] . A more advanced model, named Quantum Molecular Dynamics (QMD), was used 
to identify the various pasta structures and to determine the possible transitions between all those 
structures 

In CMD, the nucleons are represented by points in classical phase space. The QMD model, on 
the other hand, uses localised states and mimics the Pauli principle through a potential. Nevertheless, 
QMD can still be interpreted as CMD with the classical Hamilton equations of motion : 

OH . dH 

When studying phenomena at length scales smaller than the de Broglie wavelength of the system's 
constituents, CMD breaks down. Under those conditions, quantum effects stemming from the wave 
character of the nucleons predominate, and statistical many-body correlations related to the system's 
fermionic content play a fundamental role. In QMD, the effect of fermionic correlations is mimicked by 
means of a phenomenological Pauli potential. Although this method is quite successful in calculating 
energies, it fails at reproducing the single-particle occupations in momentum space for the free Fermi 
gas [HI [15]. Furthermore, in QMD, the spin and the width of the localised states are assumed to be 
time independent. 

The shortcomings of QMD for fermion systems are remedied in Antisymmetric Molecular Dynamics 
(AMD) [16J or the more general Fermionic Molecular Dynamics (FMD) [17J. In these techniques, the 
many-body state is introduced as an antisymmetrised product, i.e. Slater determinant, of localised 
single-particle states denoted as: 

|$> =i|0i(zi))g> •••(£) \M*a)). (2) 
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Here, z p represents the complex time-dependent parametrisation of the single-particle state \(j) p ). The 
time evolution of these parameters is determined by the time-dependent variational principle [18] which 
leads to the generalised Hamilton equations of motion : 

In this equation, C is the metric of the parameter space of |$) and T~L = is the Hamil- 

tonian expectation value. The matrix o represents the inverse of the overlap matrix n with n pq = (cf) p \(j) q ). 

The time-evolved many-body state is used to study the properties of the Fermi system by evaluating 
various TV-body operators. The expectation values of one- and two-body operators are calculated as 

A A 

Bi = ^(0 p |5 7 |0 g )o gp , B H = - ^2 ( ( frp ( i > r\Bii\<j) q <j)s)(o q pO sr - o qr o sp ). (4) 

pq=l pqrs=l 

Because of the determinant structure of |<&), the equations relevant to FMD can be written in matrix 
notations using the inverse overlap matrix o. Thereby, the matrices n and o play a fundamental role 
as they carry all information about the fermion statistics of the system under study. As their eigenval- 
ues can cover many orders of magnitudes, it is paramount to calculate these matrices as accurately as 
possible, i.e. analytically. 

Both AMD and FMD enjoy successes in the study of heavy-ion collisions, ground states of nuclei and 
astrophysical reaction calculations [16j HHJ [20] . However, due to the long-range character of the Pauli 
correlations, evaluating bulk fermion matter becomes a tedious task. The dimension of the matrices 
usually impedes simulations of a large number of particles, however desirable these may be for bulk 
fermion systems such as crust al matter. 



3 Bulk fermionic molecular dynamics 

When investigating the properties of bulk matter by means of large simulation volumes, the evaluation 
of expectation values becomes cumbersome. Moreover, surface effects may influence the results when a 
large fraction of the constituents lies on the surface of the simulation volume. A time-honoured method 
is to introduce a periodic structure in the simulation. When studying bulk matter using a periodic 
structure, it is crucial to make sure that the studied properties of the small but infinitely repeated pe- 
riodic system and the macroscopic system which it represents, are the same. As long as the correlation 
volume of the interactions does not exceed the simulation volume, the imposed periodicity works fine. 
However, serious problems arise in the presence of long-range correlations as e.g. induced by Coulomb 
interactions or by the Pauli exclusion principle. 

When studying bulk fermion systems in a mean-field approach using Slater determinants, the single- 
particle states are generally delocalised single-particle states fulfilling certain boundary conditions. 
These could be "periodic boundary conditions" , "Bloch conditions" or "twist- averaged boundary condi- 
tions" [211 [22]. Spatial localisations are, however, common in the neutron-star crust. To describe these 
with delocalised states, a large configuration space is required. The localised states, as used in FMD, 
would circumvent this problem but are difficult to use for bulk inhomogeneous matter. 

As the Pauli exclusion principle introduces long-range many-body correlations, it poses a major 
challenge to the study of bulk fermion systems with non-orthogonal localised states. It is indispensable 
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to fully antisymmetrise the many-body wave function representing the system. To simulate bulk matter 
with localised states, we propose to make the spatial positioning of the single-particle states periodic. 
A number of particles are placed in unit cells on a lattice 23, that is tessellating space perfectly. Each 
cell, containing A particles, can be identified by a lattice vector R = n\a\ + n 2 a 2 + n 3 a 3 with integer 
rij. The trial state of the system is thus, with T(R) the translation operator over ii, 

|$oo) = A (g) f(R) {|0!> ® • • • <g> \cf) A )} . (5) 

For |$oo), the overlap matrix n and its inverse o have infinite dimensions. Nonetheless, they exhibit 
a peculiar nested block- Toeplitz structure which can be exploited [23j . As a result of the translational 
invariance of the system, each Ax A block of the overlap matrix can be identified unambiguously with 
the lattice vector R that connects the two cells of the bra and ket states. The blocks can be evaluated 
as n pq ,R = {(j) p \T(—R)\(j) q ). Although n and o are of infinite dimension, it is possible to obtain the 
inverse overlap matrix through the following scheme 

Af(k) = V n R e- ikR , 0(k)=Af(k)-\ o R = -L f 0{k) e ikR d 3 k, (6) 

R& VBZ JBZ 

where BZ represents the first Brillouin zone of the lattice OS and k is a vector in this volume [23] . 

The strength of the proposed formalism reveals itself upon evaluating the operators in reciprocal 
space. Normally, the expectation values of the infinite fermion system would be calculated by means of 
Eqs. Q resulting in infinite sums over the block structure. These sums, however, translate into integrals 
over the first Brillouin zone which are computationally straightforward to evaluate. The expectation 
value of a one-body operator per unit-cell volume can then be computed as 

B P j = ^- / E ®iAk)OM d 3 k, B Im {k) = £ (0 p \Bjf(R)\0 q ) e ikR . (7) 

BZ Jbz pq=l Re% 

Here we assume that the operator Bj commutes with the translation operator T(R). The two-body 
operator and the metric have analogue structures [23J. 

From a computational perspective, a proper choice of the localised single-particle states helps 
speeding-up the computations. Gaussian wave packets come with the interesting feature that most 
matrix-elements can be evaluated analytically. Furthermore, the computational order for the calcu- 
lation of expectation values is the same for bulk fermion systems as for finite-sized systems (A^ 4 for 
two-body interactions), even though the former represent a system with an infinite number of particles. 
On the other hand, the equations for the bulk systems involve an extra integration over the Brillouin 
zone of the imposed lattice, increasing the computational effort. However, it has been shown that 
integrals over the Brillouin zone can be performed using only a limited number of specific fc-points [24] . 



4 Results 

In the following, we will show that the proposed technique reproduces the features intrinsic to the 
fermionic behaviour of the system. The unnormalised single-particle states are Gaussian wave packets 
of the form (x\ab) = exp{ — (x — b) 2 /(2a)} where the complex vector b represents the mean position in 
phase space and a is a complex parameter connected with the width of the wave packet. In Ref. [T7] it 
was shown that for a hundred periodically positioned wave packets, Eqs. Q reproduce the momentum 
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and spatial densities intrinsic to one-dimensional Fermi systems. In Fig. [TJ we show that our method 
also reproduces this result by placing a single particle in a unit cell of size £. It is clearly visible that 
with increasing \fajt and thus with increasing overlap of neighbouring states, hence increasing influence 
of the antisymmetrisation, the fermionic behaviour becomes apparent. For yfajt = 1, the momentum 
distribution does not show a sharp cutoff near the Fermi momentum because the periodic set of Gaus- 
sian wave packets does not constitute a complete set. For y/a/£ — >■ oo, the antisymmetry operator 
projects the Gaussian wave packets on plane waves and a sharp cutoff is reached. Even though this 




scaled position x/i scaled momentum \k\/k F 

Figure 1: Spatial density (panels (a) and (c)) and momentum density (panels (b) and (d)) of equidistant 
Gaussian wave packets with one-dimensional periodicity. One Gaussian with variance a per unit cell of 
size £. For small overlap (top) the system behaves like one of distinguishable particles while for large 
overlap (bottom) the fermionic behaviour becomes manifest. For comparison, the spatial distribution of 
a uniform Fermi gas (y/a/£ — >• oo) is presented in panel (c) (dashed line) and the momentum distribution 
of distinguishable particles is presented in panel (d) (dotted line). 

one-dimensional example clearly hints at free-fermion behaviour, it has to be assured that fermion-like 
properties of the simulations are not an artifact of the symmetry imposed by the lattice. In fact, it 
can be shown that in more dimensions the momentum distribution evolves towards the first Brillouin 
zone of the lattice [23J. This effect can be expected as, with increasing overlap, the wave packets seem 
uniform to the unit cell and the system behaves as plane waves in a periodic system. 

As we introduced periodic boundary conditions to eliminate surface effects, a successful investigation 
of bulk matter obviously requires that influences of the geometry of the chosen boundary conditions 
are negligible. In Fig. [2] we show that this is achieved for a situation where the 25 single-particle states 
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are randomly distributed in the cell. The figure shows that the spatial and momentum distribution of 
the system reflects that of a free fermion system. The perturbations on the Fermi-sphere in momentum 
space, as well as those on the uniform distribution in coordinate space can be understood as a result of 
the local clustering. Hence, a simulation with the proposed periodic structure reproduces the intrinsic 
bulk fermion behaviour, independent of the imposed periodic structure. As was indicated earlier, this 
result was hard to reach using QMD. For comparison, the bottom panels of Fig. [2] display the densities 
for distinguishable particles occupying the same wave packets. Both the spatial and momentum distri- 
butions are quite different from the indistinguishable fermion case. It clearly shows that antisymmetry 
induces Fermi motion and redistributes density in coordinate and momentum space. This long-range 
many-body correlation should not be ignored. 



Coordinate density p x i 2 /N Momentum density p k kj/N 




Figure 2: Spatial and momentum distribution of TV = 25 Gaussian wave packets with periodic boundary 
conditions. The wave packets are randomly placed in the unit cell of a square lattice with unit length 
£ and have a variance given by ^/a = 0.2£. The mean momentum of the individual wave packets is 
set to zero. The densities are normalised and coordinates are expressed in units of £ or ki = VNtt/£. 
The upper two panels depict the antisymmetrised case. The lower two panels represent distinguishable 
particles. The black circles show the positions of the centroids of the wave packets. The white circle 
represents the Fermi "sphere" of a system of free particles with the same mean density, hp = V^ttN /£. 



5 Conclusion 

In summary, we have implemented periodic boundary conditions in the FMD formalism. This allows 
one to study infinitely extended non-uniform fermionic matter with full antisymmetrisation. The intro- 
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duced spatial periodicity in the trial- state leads to a Toeplitz structure for the FMD matrices, a feature 
that can be exploited to reduce the numerical cost of the simulation. The structure of the resulting 
equations resemble those of finite FMD done on a single unit cell, however with additional integrations 
in the reciprocal lattice space. This new formalism offers perspectives to compute the expectation values 
of various operators for bulk fermion matter in a computationally feasible fashion. Although the equa- 
tions only address a finite number of particles N in a unit cell, they keep track of the Fermi statistics 
of the infinite system. We have evaluated the validity of this description with a study of the spatial 
and momentum distribution of various lattice systems and showed that our simulations convincingly 
reproduce free fermion behaviour. This shows that full antisymmetrisation can be imposed in the bulk 
FMD formalism. 

This work was supported by the Fund for Scientific Research Flanders (FWO), the Research Council 
of Ghent University and the Helmholtz International Center for FAIR (HIC for FAIR) . 
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